#
#
# LSTM weather prediction demo
# Written by: Dan R 2020
#
#

#
# Core Keras libraries
#
# import seq2seq
# from seq2seq.models import SimpleSeq2Seq
import math

from tensorflow import keras
from tensorflow.keras import Input, Model
from tensorflow.keras.layers import Embedding, Attention
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense
from tensorflow.keras.layers import LSTM, RepeatVector, TimeDistributed
from tensorflow.keras.layers import Bidirectional

# from tensorflow.contrib import seq2seq
# import seq2seq
# from seq2seq.models import SimpleSeq2Seq

#
# For data conditioning
#
from scipy.ndimage import gaussian_filter1d
from scipy.signal import medfilt

#
# Make results reproducible
#
from numpy.random import seed

seed(1)
# from tensorflow import set_random_seed
import tensorflow as tf
import pandas as pd  # 数据分析包
import numpy as np  # 提供多维数组对象的库

tf.random.set_seed(1)

#
# Other essential libraries
#
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score
from numpy import array

# Make our plot a bit formal
font = {'family': 'Arial',
        'weight': 'normal',
        'size': 10}
plt.rc('font', **font)

#
# Set input number of timestamps and training days
#
n_timestamp = 5
train_days = 1200  # number of days to train from
testing_days = 400  # number of days to be predicted
n_epochs = 20
filter_on = 1
train_model = 0
# x_or_y = 0     # 1x 0y
user_id = '001'  # 1 3 5 8 11 14 15 38

x_delta = np.zeros((n_timestamp, testing_days - n_timestamp - 1))
y_delta = np.zeros((n_timestamp, testing_days - n_timestamp - 1))

# seq2seq
n_in = 5
n_out = 5
encoded_length = 1
batch_size = 7


# 创建一个差分序列
def difference(dataset, interval=1):
    diff = list()
    # print("len(dataset)",len(dataset))
    for i in range(interval, len(dataset)):
        value = dataset[i] - dataset[i - interval]
        diff.append(value)
    return diff


# 差分的逆转换
def inverse_difference(last_ob, value):
    return value + last_ob


#
# Select model type
# 1: Single cell
# 2: Stacked
# 3: Bidirectional

#
model_type = 5

url = "C:/Users/admin/Desktop/WeatherLSTM-master/NewYork/"
dataset = pd.read_csv(url + 'NewYork_30sec_' + user_id + '.csv', error_bad_lines=False)
print(dataset)
if filter_on == 1:
    dataset['X'] = medfilt(dataset['X'], 3)
    dataset['X'] = gaussian_filter1d(dataset['X'], 1.2)
    dataset['Y'] = medfilt(dataset['Y'], 3)
    dataset['Y'] = gaussian_filter1d(dataset['Y'], 1.2)
# print("dataset[X]=",dataset['X'])
# print("dataset[Y]=",dataset['Y'])
#
# Set number of training and testing data
#
train_set = dataset[0:train_days].reset_index(drop=True)
test_set = dataset[train_days: train_days + testing_days].reset_index(drop=True)
for x_or_y in range(2):
    if x_or_y:
        training_set = train_set.iloc[:, 1:2].values
        testing_set = test_set.iloc[:, 1:2].values
        diff_train = difference(training_set)
        diff_test = difference(testing_set)
        # print("diff_train", len(diff_train))
        # print("diff_test", len(diff_test))
    else:
        training_set = train_set.iloc[:, 2:3].values
        testing_set = test_set.iloc[:, 2:3].values
        diff_train = difference(training_set)
        diff_test = difference(testing_set)
        # print("diff_train_y", len(diff_train))
        # print("diff_test_y", len(diff_test))
    #
    #
    # Normalize data first
    #
    sc = MinMaxScaler(feature_range=(0, 1))

    training_set_scaled = sc.fit_transform(diff_train)
    testing_set_scaled = sc.fit_transform(diff_test)


    #
    # Split data into n_timestamp
    #
    def data_split(sequence, n_timestamp):
        X = []
        y = []
        for i in range(len(sequence)):
            end_ix = i + n_timestamp
            if end_ix > len(sequence) - 1:
                break
            # i to end_ix as input
            # end_ix as target output
            seq_x, seq_y = sequence[i:end_ix], sequence[end_ix]
            X.append(seq_x)
            y.append(seq_y)
        return array(X), array(y)


    X_train, y_train = data_split(training_set_scaled, n_timestamp)
    X_train = X_train.reshape(X_train.shape[0], X_train.shape[1], 1)
    X_test, y_test = data_split(testing_set_scaled, n_timestamp)
    read_x_test = X_test
    # print("len X_test", len(X_test))
    X_test = X_test.reshape(X_test.shape[0], X_test.shape[1], 1)
    # print("X_train=",X_train)
    # print("X_train.shape[0]=",X_train.shape[0]) 1195
    # print("X_train.shape[1]=",X_train.shape[1]) 5

    # X_train = X_train.reshape(X_train.shape[0], X_train.shape[1], 2)

    # X_test = X_test.reshape(X_test.shape[0], X_test.shape[1], 2)
    if model_type == 5:
        model = Sequential()
        model.add(LSTM(50, input_shape=(X_train.shape[1], 1)))
        model.add(RepeatVector(n_out))
        model.add(LSTM(50, return_sequences=True))
        model.add(TimeDistributed(Dense(1)))

    #
    # Start training
    #
    if train_model:
        model.compile(optimizer='adam', loss='mean_squared_error')
        model.summary()
        history = model.fit(X_train, y_train, epochs=n_epochs, batch_size=32)
        loss = history.history['loss']
        epochs = range(len(loss))
        if x_or_y:
            model.save("./model/" + "seq2seq_model_" + user_id + "_x.h5")
        else:
            model.save("./model/" + "seq2seq_model_" + user_id + "_y.h5")
    else:
        if x_or_y:
            model = keras.models.load_model("./model/" + "seq2seq_model_" + user_id + "_x.h5")
        else:
            model = keras.models.load_model("./model/" + "seq2seq_model_" + user_id + "_y.h5")
    #
    # Get predicted data
    #
    y_predicted = model.predict(X_test)
    # print("y_predicted=", y_predicted)
    y_predicted = np.reshape(y_predicted, (-1, 1))
    # print("长度y_predicted=", y_predicted.shape[0])
    b = np.zeros((5, 394))
    for a in range(394 * 5):
        if a % 5 == 0:
            b[0][int(a / 5)] = y_predicted[a]
        elif a % 5 == 1:
            b[1][int(a / 5)] = y_predicted[a]
        elif a % 5 == 2:
            b[2][int(a / 5)] = y_predicted[a]
        elif a % 5 == 3:
            b[3][int(a / 5)] = y_predicted[a]
        elif a % 5 == 4:
            b[4][int(a / 5)] = y_predicted[a]
    # print("可分离出来每一步的b=", b)
    #
    # 'De-normalize' the data
    #
    y_predicted_descaled = sc.inverse_transform(b)
    y_train_descaled = sc.inverse_transform(y_train)
    y_test_descaled = sc.inverse_transform(y_test)
    # print("len y_predicted_descaled",len(y_predicted_descaled[0]))
    # print("len y_test_descaled",len(y_test_descaled))
    # print("y_predicted_descaled",y_predicted_descaled)
    # print("y_test_descaled",y_test_descaled)
    for i in range(len(y_predicted_descaled[0])):  # 反差分
        y_test_descaled[i] = y_test_descaled[i] + testing_set[n_timestamp + i]
        for a in range(5):
            y_predicted_descaled[a][i] = y_predicted_descaled[a][i] + testing_set[n_timestamp + i]
    # print("y_predicted_descaled", y_predicted_descaled)
    # print("y_test_descaled", y_test_descaled)
    # print("y_predicted_descaled-y_test_descaled", y_predicted_descaled - y_test_descaled)
    # print("pre-true",y_predicted_descaled-y_test_descaled)
    y_pred = y_predicted.ravel()
    y_pred = [round(yx, 2) for yx in y_pred]
    y_tested = y_test.ravel()
    for i in range(testing_days - n_timestamp - 1):  # 对于与上一时刻一样的点，直接赋原值
        if i >= 2:
            if y_test_descaled[i - 2] == y_test_descaled[i - 1]:
                for a in range(5):
                    y_predicted_descaled[a][i] = y_test_descaled[i - 1]
    predict_sum = 0
    for i in range(testing_days - n_timestamp - 1):
        predict_sum = predict_sum + abs(y_predicted_descaled[0][i] - y_test_descaled[i])
        if x_or_y:
            for a in range(5):
                x_delta[a][i] = abs(y_predicted_descaled[a][i] - y_test_descaled[i])
        else:
            for a in range(5):
                y_delta[a][i] = abs(y_predicted_descaled[a][i] - y_test_descaled[i])
    predict_delta_advange = predict_sum / (testing_days - n_timestamp - 1)
    if x_or_y:
        print("predict_delta_advange_x", predict_delta_advange)
        y_test_descaled_x = y_test_descaled
        y_predicted_descaled_x = y_predicted_descaled
    else:
        print("predict_delta_advange_y", predict_delta_advange)
        y_test_descaled_y = y_test_descaled
        y_predicted_descaled_y = y_predicted_descaled
    # print("y_test_descaled", y_test_descaled)
    # print("len y_test_descaled", len(y_test_descaled))
predict_sum_xy = 0
for i in range(testing_days - n_timestamp - 1):
    predict_sum_xy = predict_sum_xy + np.sqrt(x_delta[0][i] ** 2 + y_delta[0][i] ** 2)
predict_delta_advange_xy = predict_sum_xy / (testing_days - n_timestamp - 1)
print("predict_delta_advange_xy", predict_delta_advange_xy)

'''
#
# Show results
#
plt.figure(figsize=(8, 7))

plt.subplot(3, 1, 1)
plt.plot(dataset['X'], dataset['Y'], color='black', linewidth=1, label='True value')
plt.ylabel("Y")
plt.xlabel("X")
plt.title("All data")

# plt.subplot(3, 2, 3)
# plt.plot(y_test_descaled, color='black', linewidth=1, label='True value')
# plt.plot(y_predicted_descaled, color='red', linewidth=1, label='Predicted')
# plt.legend(frameon=False)
# plt.ylabel("X")
# plt.xlabel("Time")
# plt.title("Predicted data (n Times)")

plt.subplot(3, 2, 3)
plt.plot(y_test_descaled_x, y_test_descaled_y, color='black', linewidth=1, label='True value')
plt.plot(y_predicted_descaled_x[0], y_predicted_descaled_y[0], color='red', linewidth=1, label='Predicted')
plt.legend(frameon=False)
plt.ylabel("X")
plt.xlabel("Y")
plt.title("Predicted location")

plt.subplot(3, 2, 4)
plt.plot(y_test_descaled_x, color='black', linewidth=1, label='True value')
plt.plot(y_predicted_descaled_x[0], color='red', label='Predicted')
plt.legend(frameon=False)
plt.ylabel("X")
plt.xlabel("Time")
plt.title("Predicted data (first 400 Times)")
# plt.subplot(3, 2, 4)
# plt.plot(y_test_descaled_x[0:75],y_test_descaled_y[0:75], color='black', linewidth=1, label='True value')
# plt.plot(y_predicted_descaled_x[0:75],y_predicted_descaled_y[0:75], color='red', label='Predicted')
# plt.legend(frameon=False)
# plt.ylabel("X")
# plt.xlabel("Time")
# plt.title("Predicted data (first 75 Times)")

# plt.subplot(3, 3, 7)
# plt.plot(dataset['Y'], color='black', linewidth=1, label='True value')
# plt.ylabel("Y")
# plt.xlabel("Time")
# plt.title("All data")
# plt.subplot(3, 3, 7)
# plt.plot(epochs, loss, color='black')
# plt.ylabel("Loss (MSE)")
# plt.xlabel("Epoch")
# plt.title("Training curve")
plt.subplot(3, 3, 7)
plt.plot(y_test_descaled_y, color='black', linewidth=1, label='True value')
plt.plot(y_predicted_descaled_y[0], color='red', label='Predicted')
plt.legend(frameon=False)
plt.ylabel("Y")
plt.xlabel("Time")
plt.title("Predicted data (first 400 Times)")

plt.subplot(3, 3, 8)
plt.plot(y_test_descaled_x - y_predicted_descaled_x[0], color='black')
plt.plot(y_test_descaled_y - y_predicted_descaled_y[0], color='red')
plt.ylabel("Residual")
plt.xlabel("Time")
plt.title("Residual plot")

# plt.subplot(3, 3, 9)
# plt.scatter(y_predicted_descaled, y_test_descaled, s=2, color='black')
# plt.ylabel("Y true")
# plt.xlabel("Y predicted")
# plt.title("Scatter plot")

plt.subplots_adjust(hspace=0.5, wspace=0.3)
plt.show()

mse = mean_squared_error(y_test_descaled, y_predicted_descaled[0])
r2 = r2_score(y_test_descaled, y_predicted_descaled[0])
print("mse=" + str(round(mse, 2)))
print("r2=" + str(round(r2, 2)))
'''
